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We study the spectrum of a random matrix, whose elements depend on the Euclidean distance 
between points randomly distributed in space. This problem is widely studied in the context of the 
Instantaneous Normal Modes of fluids and is particularly relevant at the glass transition. We intro- 
duce a systematic study of this problem through its representation by a field theory. In this way we 
can easily construct a high density expansion, which can be resummed producing an approximation 
to the spectrum similar to the Coherent Potential Approximation for disordered systems. 

C/2 ' I. INTRODUCTION 

The theory of random matrices has found applications in many branches of physics [h| . The most 
developed theory concerns matrices where the matrix elements are either independent random vari- 
ables, as in Gaussian ensembles, or are taken with a statistical distribution which is invariant under 
some symmetry group. However in many physical applications, from vibration spectra of glasses 
plpl to instantaneous normal modes in liquids H] , electron hopping in amorphous semiconductors H 
or combinatorial optimization problems M, one needs to compute the spectrum and the eigenstate 
properties of some random matrices which are of a different type: The disorder is due to the random 
positions of N points, and the matrix elements are given by a deterministic function of the distances 
between the points. We shall call these matrices Euclidean Random Matrices(ERM). 

Specifically we want to study the following mathematical problems: Consider N points in a volume 
V of a d-dimensional Euclidean space. For a given sample, characterized by the positions Xi of the 
iV points (i g {1...JV}), we want to study the properties of the N x N random matrices M defined 
as 
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where it is a real parameter which enables us to interpolate between the two most interesting cases 
u — 0,1. The case where u = is the simplest mathematical problem with Euclidean-correlated 
matrix elements, the case where u — 1 is the natural problem which appears when studying for 
instance vibration modes of an amorphous solid, instantaneous normal modes of a liquid, or random 
master equations. The main difference is that when u = 1 the matrix M has fluctuating diagonal 
terms, tailored in such a way that JJ . My = 0: the vector with all components equal to one is an 
eigenvector with zero eigenvalue, which expresses global translation invariance. 

To fully specify the problem one needs to characterize the probability distribution of the random 
points Xi, as well as the function f(x). In this paper we shall concentrate on the case where the 
points are uniformly distributed in a cubic box of size L — V 1 ' , without correlations. In many 
applications one will have to generalize the problem to the case where the x^s are correlated, 
including a short distance repulsion. This is a much more complicated problem which we shall not 
address here. As for the function f(x), we shall assume that it depends only on the distance \x\ 
and that it decays fast enough at large argument. In particular we shall assume that the Fourier 
transform f(k) = J dxe 1 x fix) is a well defined function at all k. 

The general problem of ERM theory is to understand the statistical properties of the eigenvalues 
and the corresponding eigenvectors of M in the large N limit (taken at fixed density p — N/V). 
This is certainly a very rich problem. In particular one can expect that the eigenvectors will exhibit 
in dimension d > 3 some localized and delocalized regimes separated by a mobility edge H|. Here 
we wish to keep to the much simpler question of the computation of the spectrum of M . We shall 
develop a field theory for this problem, check it at high and low densities, and use a Hartree type 
method to approximate the spectrum in the u — case. 



II. FIELD THEORY 

The spectrum can be computed from the resolvent: 



where the overline denotes the average over the positions Xi- It is possible to write down a field 
theory using a replica approach. We shall compute Hjv = det(z — M)~ n / 2 , and deduce from it the 
resolvent by using the replica limit n — » 0. It is easy to show that one can write Hjv as a partition 
function over replicated fields </>", where i g {1...JV}, a g {l...n}: 



n^/nn-w 



dxi 

2 = 1 a=l 

ex p ( -fX^) 2+ 2 5Z ^^ _ x ^^^ _ l^w) 2 ^/^ - ^ - ) ) ■ ( 3 ) 

Let us introduce the bosonic fields ip a {x) = ^ i=1 0"<5(a; — £;) and x( 2; ) — Xy- (fit) $(x — Xi), 
together with their respective Lagrange multiplier fields ip a (x) and x(x). One can integrate out the 
<j> variables, leading to the following field theory for Hjv: 
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It is convenient to go to a grand canonical formulation of the disorder: we consider an ensemble 
of samples with varying number of points, and compute the grand canonical partition function 
Z(a) = X^iv = o ^ncc n /N\ which is equal to: 
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Z= / L>[Va,Va,X.x]exp(5o + S'i) ; S 1 =aA. (6) 

Notice that we can also integrate out the ip field thus replacing So by So, where 

S'o = ~tr log / + - / dx x (x)x(x) + -^2 / dxd V ^ i x ) f' 1 {x ~ y)4> a (y) , (7) 
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where log / is the logarithm of / considered as an integral operator and / _1 is the integral operator 
which is the inverse of /. 

The expression (pi) is our basic field theory representation. We shall denote by brackets the 
expectation value of any observable with the action So + Si. As usual with the replica method 
we have traded the disorder for an interacting replicated system. The basic properties of the field 
theory are related to the properties of the original problem in a straightforward way. The average 
number of particles is related to a through N — a{A), so that one gets a — N in the n — > limit. 
From the generalized partition function Z, one can get the resolvent R(z) through: 

i?(z) = _l im A^ (8) 
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III. HIGH DENSITY EXPANSION 

Let us first show how this field theory can be used to derive a high density expansion. We shall 
rescale z as z = pz, and rescale the \ fields as: x( x ) = pc(x), x( x ) — C ( x )/P- The interacting part 
of the action, Si, can then be expanded as: 
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where we have introduced: 
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K X ,V)= / drf(x,r)p(r)f(r,y) 



(9) 
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Performing the quadratic c functional integral, one finds (here and in the following we drop all the 
irrelevant -z independent- constants): 
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dxp(x)~- / dxdy e(x)h x (x, y)e(y) - - ^ / dxdy ip a (x)G 1 (x,y)i> a (y) 
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where G J (a;, j/) = <5(x — j/)/(£ + c(o;)) — / 1 (x,y) is the propagator of the Tp field, and the field e(x) 
is defined as: 

e(aO = — f c(x) - u / dy/(x - j/)/ti(y) j = — f c(x) - u/(0) - — / dy/(x - y) log(z + c(y)) 

We then change variables from c to e, which involves a Jacobian which is easily computed to order 
n. The ip integral is Gaussian, and to leading order in p, the e field can be neglected in G, as well 
as in h. The only term in which e is relevant turns out to be J dxp(x), which is simplified by an 
expansion to order e 2 . Gathering the various contributions one gets 
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where the first two terms come from the expansion of J dxp(x), the third one is the Jacobian, and 
the fourth one is the contribution of the determinant from the tp fluctuations. 
This gives for the resolvent: 
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This result is valid in the large p limit, whenever the resolvent parameter z is scaled as z — pz. It 
can be checked directly. We consider the original form (g) of the resolvent and expand it in a 1/z 
series. At order l/z p+1 , the first two leading terms in p are: 
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where Cf } = F(p + l)/(F(p' + l)T(p — p' + 1)). One can see that in the large p limit the series 
simplifies if we first scale z proportionally to p. Writing z = pz, we can resum the leading terms in 
the series, which gives back the same resolvent as (Ml). 



One can study with this method the eigenvalue density for eigenvalues |A| ~ O(p), by taking 
z — A + ir\ and computing the imaginary part of the resolvent in the small r\ limit. For p — » oo one 
would get the trivial result for the eigenvalue densityP(A) = S(X + puf(0)). Including the leading 
large p correction which we have just computed, we find that P(X) develops, away from the peak at 
A ~ — pu/(0), a component of the form: 

P(X) ~ - p J(dk)S(X - p(f(k) - u/(0))) . (16) 

The result (hq) is the one that one could guess using the following simple argument: introducing 
4>j = exp(ikxj), one gets: 

^Miitj = Xifc , Ai = 5^e _ ** ( "«-"' ) /(a!i-a! i )-tt5^/(a!i-« m ). (17) 
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In the large density regime it is reasonable to use the central limit theorem to approximate the 
above sum, which would show that cf> is near to an eigenvector, with eigenvalue p(f(k) — u/(0)). 
Notice that for this argument to hold, the discrete sum giving Aj in fll7| ) must sample correctly the 
continuous integral. This will be the case only when the density p is large enough that the phase 
—ik(xi — Xj) doesn't oscillate too much from one point Xj to a neighbouring one. This imposes that 
the spatial frequency |fc| be small enough: \k\ <C p 1 ' 4 . This same condition is present in the field 
theory derivation. We assume that f(k) decreases at large k, and we call Hm the typical r ang e of 
k below which f(k) can be neglected. Let us consider the corrections of order p~ 2 in eq. (|l4[). It 
is clear that, provided z is away from — w/(0), the ratio of the correction term to the leading one is 
of order kfjp^ 1 , and the condition that the correction be small is just identical to the previous one. 
The large density corrections near to the peak A ~ —puf(0) cannot be studied with this method. 



IV. LOW DENSITY EXPANSION 

The low density expansion is also easily performed from the field theoretic representation. Since 
the interaction a.A is proportional to p in the n — > limit, we can expand in powers of a: 

Z = j D[^ a , X ,x\e s 'o 
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At first order we can perform the \ integral, which fixes x(x) = uf(x — xo) and we get: 

Z = l + p(^-) (detA-i)-™/ 2 (19) 

where K\ is the operator: K±(x,y) — S(x — y) — S(x — xa)5(y ~ xo)/(0)/(z + u/(0)). Using (g), 
one finds R(z) = l/(z — (1 — u)/(0)), which is obviously the leading result at very low densities, 
such that the points are isolated. At second order, the expansion (n8|) involves a double integral 
over points xq,xi. The x integral fixes x( x ) — uf(x — Xo) + uf(x — xi). One needs to study the 
determinant of the operator K2 defined as: 

K 2 = 5(x -y)- /( f 7 V \ (S(x - xo) + 5(x - xi)) (S(y - x ) + 5(y - xi)) . (20) 

z + X\x) 

This gives after a simple computation the order p correction to the resolvent: 

2 J dT { z - [(1 - u)f(0) + (1 - u)f(r)] + z - [(1 - u)f(0) - (1 + u)f(r)] ~ z - (1 - «)/(0) J ' (21) 

This corrects the leading small p formula by adding the contribution due to pairs of points, a distance 
r apart, which are isolated from all other points. This expansion can be easily carried out to higher 
order, the order p k involving the computation of a kd dimensional integral. 



V. VARIATIONAL APPROXIMATION 

In order to try to elaborate a general approximation for the spectrum, which interpolates between 
the high and low density limits, we have used a standard Gaussian variational approximation in the 
field theory representation, which appears under various names in the literature, like the Random 
Phase Approximation. We have developed it only in the case u = 0. The field theory then simplifies 
since the \ an d X fields can be integrated out. Changing %f) — > iip, we obtain: 

L>[^ a ]exp(S* u=0 ) , (22) 

where 
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which approximates the full interacting problem, using the fact that the variational free energy 
F v =< 5*11=0 > v —logZ v should be minimized. (Here Z v = J D[tp] exp(SV,), and the expectation 
values < . >„ are meant as Boltzmann like averages with the measure exp(5„)). The variational 
free energy is easily obtained: 
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where the operator K is defined as: 
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where H a b — G a t(r = 0) = J(dk)G a b(k) is a n x n matrix, and the trace Tr n is a trace in replica 
space. This gives for the variational free energy: 
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variational approximation to the spectrum for any values of / and the density. In sect. VI we shall 
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We look for the best quadratic action 

S v = -(1/2) J2 I dxdyG- b \x, y)i> a (x)^ b (y) (24) 
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K ab( x ,y) =Gj{x,y)- -5 ab S(x ~ x )S(y - x ) ■ (26) 

In (Eq), both G and K are considered as operators in x space and in replica space. Keeping 
to translation invariant variational actions, the result is easily expressed in terms of the Fourier 
transform G a t(k) of the function G a b(r). One finds that: 
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§ = \J{dk) Tr J^ k) - M exp (-lrr„log(l - H/z)) - ~ j (dk)Tr n logG(fc) (28) 



We have found a solution of the stationarity equations, dF v /dG a b(x) — 0, of the form G a b(x) — 
S a bG(x), where the Fourier transform G satisfies the self consistency equation: 

e ^-T%s- c - .. I ^em (29) 

In terms of this function, the variational free energy density per replica is: 

^ = iy (d fe)|M + £ log (l-lJ(dk)G(k)\ -\J(dk) log G(k) (30) 

and the resolvent is obtained from: 

R( Z ) = Ti 77 = = 7. = (31) 
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Formulas (E9lpd,pij) provide a closed set of equations which allow us to compute the Gaussian 



compare this approximation to some numerical estimates of the spectrum for various functions / and 
densities. It is easily checked that the variational gives the correct spectrum to leading order in the 
large p expansion (see (114)) and in the small p expansion, which is not a surprise since we had seen 
that these leading orders are given by purely Gaussian theories. Of course it fails to give the exact 
result beyond leading order, but it provides a reasonable interpolation between these extremes. 

Another partial resummation of the p expansion can also be done in the following way: if one 
neglects the triangular-like correlations between the distances of the points (an approximation which 
becomes correct in large dimensions, provided the function / is rescaled properly with dimension), 
the problem maps onto that of a diluted random matrix with independent elements. This problem 
can be studied explicitly using the methods of [p|-|l2|. It leads to integral equations which can be 
solved numerically. The equations one gets correspond to the first order in the virial expansion in 
eq. (0) , in which one introduces as variational parameter the local probability distribution of the 
field <j>. The detailed computations involving this other approximation scheme are left for future 
work. 



VI. NUMERICAL SIMULATIONS 

For a function f(x) which is positive, as well as its Fourier transform f(k) , simple bounds on 
the spectrum can be derived in the two cases u — and u — 1. In the u — case, calling ipi a 
normalised eigenvector ipi of the matrix M defined in (hi), with eigenvalue A, one has: 

^2 ^f( Xi ~ x i)^i = A ' ( 32 ) 

ij 

and the positivity of the Fourier transform of / implies that A > 0. In the u — 1 case, the eigenvalue 
equation implies that: 
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Summing over i, using the fact that Mij is positive for i 7^ j, together with the constraint ^ . Mij = 
0, one gets: 

^(|Af„|-|A-M J3 |)|^|>0 (34) 
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which implies that there exists at least one index j such that \Mjj\ — |A — Mjj\ > 0, and therefore 
A < 0. Furthermore, in this u = 1 case, there exists one eigenvector (the uniform one) with zero 
eigenvalue. 

In the u — case, we have studied numerically the problem in dimension d — 3 with the Gaussian 
function f(x) — (2ty)~ 3 ' 2 exp(— x 2 /2). In this Gaussian case the high density approximation gives 
a spectrum 

p w~^Ki to8 s) 1/a ^- A) (35) 

Notice that this spectrum is supposed to hold away from the small A peak, and in fact it is not 
normalizable at small A. 

The variational approximation computation described in the previous section is more involved. 
From ( p^ , pfj|j3l] ) , one needs to solve, given z — A — ie, the following equations giving the complex 
function of z, C(z) which we write as C — a + ib: 
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One needs to find a solution in the limit where e — > 0. We have done this as follows: For a given 
value of a, we search for the values of 6 which solve the second of eqs.(B6|) with e = 0. Depending on 
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the value of a, there may exist zero, one, or two solutions in b. For each such solution we compute 
the value of A from the first of qs.(pq), and the corresponding density of states is given by b/(pw). 
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Density of eigenvalues of a Euclidean Random Matrix in three dimensions, density p — 1. The function / 
(27r) -3 ' 2 exp(— x 2 /2), and the matrix is defined from eq.(QJ) with u — 0. The full line is the result of a numerical 
simulation with N = 800 points, averaged over 100 samples. The dashed (red) line is the result from the high density expansion. 
The dash-dotted (green) line is the result from the Gaussian variational approximation (RPA) to the field theory. 
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FIG. 2. Density of the logarithm (in base 10) of the eigenvalues of a Euclidean Random Matrix in three dimensions, density 
p — 1. The function / is fix) — (27T) -3 ' 2 exp(— x 2 /2), and the matrix is defined from eq.(nl) with u = 0. The full line is the 
result of a numerical simulation with N = 800 points, averaged over 100 samples. The dashed (red) line is the result from the 
high density expansion. The dotted (green) line is the result from the Gaussian variational approximation (RPA) to the field 
theory. 



In fig. Ql|J2[) , we plot the obtained spectrum, averaged over 100 realizations, for N = 800 points 
at density p — 1 (We checked that with N = 600 points the spectrum is similar). Also shown 



are the high density approximation (pq), and the result from the variational approximation. We 
see from fig.([jj) that the part of the spectrum A 6 [0.2, 1.5] is rather well reproduced from both 
approximations, although the variational method does a better job at matching the upper edge. 
On the other hand the probability distribution of the logarithm of the eigenvalues (figfl) makes it 
clear that the high density approximation is not valid at small eigenvalues, while the variational 
approximation gives a sensible result. One drawback of the variational approximation, though, is 
that it always produces sharp bands with a square root singularity, in contrast to the tails which 
are seen numerically. 

In figM we plot the obtained spectrum, averaged over 200 realizations, for N = 800 points at 
density p = 0.1 (We have checked that there is no substantial variation of the plot when going from 
N — 600 to N — 800). Also shown are the low density approximation (Ell) , and the result from the 
variational approximation. We see from fig.(H) that this value of p = 0.1 is more in the low density 
regime, and in particular there exists a peak around A = /(0) due to the isolated clusters containing 
small number of points. The variational approximation gives the main orders of magnitude of the 
distribution, but it is not able to reproduce the details of the spectrum, in particular the peak due to 
small clusters. On the other hand the leading low density approximation, which is not normalizable, 
gives a poor approximation at this intermediate density. 
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Density of eigenvalues of a Euclidean Random Matrix in three dimensions, density p 
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f(x) — (2tt)~ 3 ' 2 exp(— x 2 /2), and the matrix is defined from eq. (Ill) with u = 0. The full line is the result of a numerical 
simulation with N — 800 points, averaged over 200 samples. The dashed (red) line is the result from the low density expansion. 
The dash-dotted (green) line is the result from the Gaussian variational approximation (RPA) to the field theory. 



We now turn to the u — 1 case. We have run some simulation in d — 3 with the same Gaussian 
/ of width one, at density p — 1. As seen in fig.|4 the direct simulations with periodic boundary 
conditions, and sizes TV = 1000 to 1400, show roughly the same density of states, with a very broad 
peak around A = — 1 = — p/(0). However they disagree in the tail of the spectrum near to the zero 
eigenvalue. In this region the finite size effects are more pronounced because of the gap opening at 
momenta less than 2ir(N/p)~ 1 ' 3 , which is rather large for the accessible values of N. In order to 
get a more precise result, less sensitive to the finite sizes, in this region of the spectrum, we have 
run some simulations which deal with an infinite set of points. Starting from a given sample of N 
points in a box of side L — (N/p) 1 ' 3 , one imagines buiding up an infinite cubic lattice, for which 
the unit cell is the original box. Bloch's theorem tells that there exist N bands, and the eigenstates 
of the infinite system are the products of a periodic function of period L (in each direction) times 
a plane wave e lpXj . For a given value of the momentum p £ [— tt/L, ir/L] 3 , the N eigenvalues are 
those of the modified matrix: 
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(37) 



Here dkj is the distance between Xi and Xj with periodic boundary conditions, defined as the 
minimum over all translation vectors t (translations of length multiple of L in each direction) of 
\xi — Xj — t\. The phase 4>ij is a d-dimensional vector. Its component in direction fi is related to 
the optimum translation t through: (j>£- — t^/L. Given p, the problem is thus mapped to finding 
the spectrum of a L d x L d matrix M^ p > , with randomly twisted boundary conditions. The resulting 
matrix is hermitian and can be diagonalized by standard library routines. We have averaged the 
spectra over many samples, choosing for each sample one random Bloch momentum. This method 
|13] allows to access the small eigenvalue region of the spectra. As we see from figM this whole 
region is well approximated by the high density approximation (|16[), but this approximation fails to 
reproduce the broad peak of the spectrum, as expected. 
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FIG. 4. Density of eigenvalues of a Euclidean Random Matrix in three dimensions, density p — 1. The function / is 
f(x) = (2-7T-) -3 ' 2 exp(— x 2 /2), and the matrix is defined from eq.(hl) with u — 1. The right hand graph gives an enlargement 
of the spectrum at small eigenvalues. The full line is the result of a numerical simulation with quenched-twisted boundary 
conditions (see text), with N — 500 points, averaged over 500 samples. The dashed (red) line is the result from the high 
density expansion. The dash-dotted (green) line is the result of a simulation with untwisted periodic boundary conditions, with 
N = 1000 points averaged over 50 samples. The dotted (blue) line is the result of a similar simulation, but with N = 1400 
points. 



To summarize, we have introduced a family of random matrices with Euclidean correlations, and 
developed for them a field theory representation, as well as some systematic expansions giving some 
properties of the eigenvalue density. The specra can be rather well approximated in the simplest 
case u — using some variational approximation to the field theory. In the most interesting case 
u = 1, the situation is more complicated and we could get only the behaviour of the spectrum at 
small eigenvalues, using the high density expansion. It would certainly be interesting to generalize 
the variational approximation in order to treat also this u = 1 case. 
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